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Abstract 

An analytical function is derived that exactly describes the amplification process due to a series of 
discrete, Poisson-like amplifications like those in a photo multiplier tube (PMT). A numerical recipe 
is provided that implements this function as a computer program. It is shown how the program can 
be used as the core-element of a faster, simplified routine to fit PMT spectra with high efficiency. The 
functionality of the method is demonstrated by fitting both, Monte Carlo generated and measured 
PMT spectra. 



1 Introduction 



2 An Analytical Function 



In September 1999, the LHCb RICH group tested 
Hamamatsu's 64-Multi-anode Photo Multipher 
Tubes as a possible photodetector choice for the 
LHCb RICH detector jRT(^nn| . During the data 
analysis, the need for an accurate model of the out- 
put of a PMT arose that could be fitted to the 
measured pulse height spectra, mainly in order to 
have a precise estimate of the signal lost below the 
threshold cut. In order to perform a fit to the spec- 
tra, an analytical function is needed that can be 
calculated reasonably quickly by a computer. 

Such a function is derived here. First (section EJ, 
an analytical expression is derived that describes 
the output of a PMT. The starting assumption is 
that the number of photoelectrons per event, as 
well as the number of secondary electrons caused 
by each primary electron at each stage of the dyn- 
ode chain, are well described by Poisson distribu- 
tions. Furthermore it is shown how this expres- 
sion can be adapted to avoid some of the numerical 
problems arising in the original expression, so that 
it can be calculated by a computer. A complete 
numerical recipe is given and a FORTRAN imple- 
mentation of the program is listed in appendix lAl 
This expression can of course be used to calculate 
any "snowball" like effect described by a series of 
Poisson distributions. 

In section 13 it is described how the exact expres- 
sion derived in the first part can be used as the cen- 
tral element of a faster, approximate function, and 
how the number of parameters can be reduced mak- 
ing reasonable assumptions, so that fitting a large 
number of spectra in a finite time becomes feasible. 
This is then adapted to describe the digitised out- 
put of laboratory read-out electronics, rather than 
the number of electrons at the end of a dynode 
chain. 



2.1 The Electron Probability Distri- 
bution 

In the following, an expression for the number of 
photoelectrons at the end of a dynode chain of a 
PMT is derived. The number of incident photons, 
and hence of photoelectrons produced in the cath- 
ode, is assumed to follow a Poisson distribution. 
This is appropriate for the testbeam data where 
PMTs were used to detect Cherenkov photons gen- 
erated by a particle traversing a dielectric. With 
a mean number of photoelectrons produced in the 
cathode of Ai , the probability to find ki photoelec- 
trons arriving at the first dynode is: 



A.A^-^ 



P(M = e-^ 



(1) 



The probability to find ^2 electrons after the first 
dynode is the sum over all values for fci of the prob- 
abilities P(fci), each multiplied by the probability 
that the dynode returns fc2 electrons given that fci 
arrive: 



P(fc2)= ^F(fci).F(fc2|fci). 



(2) 



fci=0 



Each of the fci electrons produces a Poisson- 
distributed response from the dynode with mean 
A2 where A2 is the gain at the T'' dynode; all fci 
electrons together produce a response distributed 
according to the convolution of fci Poisson distribu- 
tions, each with mean A2. This results in a single 
Poisson distribution with mean A2 • fci : 



P(fc2|fci)=e 



-Aofc 






(3) 



Hence the probability to find fc2 electrons after the 
first dynode is given by: 

P(fc2)=f:P(fci).e-^^'=^(Ml)!!. (4) 



fci=0 



hi 



Inserting the right-hand side of equation ^ for 

This approximate function is used in section El of P(fci) yields, after manipulation: 

the paper to fit Monte Carlo generated spectra as fe 00 x 

An^ — M.^- 2 



well as real data, demonstrating the validity of the 
method. 



F(fc2) = e 



- ^-^1 
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fci! 



fcr- (5) 



Generalising this for n — 1 dynodes yields: 
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(9) 



(6) Now each summation can be carried out in turn, 
starting with that over kn-i'- 
Each term in equation El is of the form of an 

k 

exponential series, i.e. ^, except for the last 
term with the summation parameter A:„_i, which 
appears as xr^'^"- This term can be expressed in 
terms of the fc„th derivative of e^'^""^ with respect 
to the new variable y at y = 0: 
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(10) 



then over kn-2'- 



d^n 



^yk„- 



dy 



fe„ 
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(7) 
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(11) 



Now the last term in equation can be written as 



(A„. 


^ — An u \fcii-l 
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and so on. After performing all these summations, 

the probability of finding fc„ electrons after n — 1 

(8) dynodes, with gains A2, . . . , A„, starting off with an 

average of Ai photo electrons arriving at the first 



Using El equation El can be re-written as: dynode, is given by: 



P(fc„) = 


exp(a;i exp(a;2 exp(x3 • • 

exp(a;„_iexp(2/)) •••))) 


y=o 


with Xi E 


= X,e-^'+\ 


(12) 



which in turn gives a recursive formula for the mth 
derivative: 



^™)^^ ™-lU.),^^.), (17) 

3=0 \ -^ / 



2.2 Calculating P(A;„) 



with /«) = /„ VjelN. 



With this expression, equation^lcan finaUy be cal- 
culated, by starting with /n(0) = 1 and calculating 
f^' subsequently for all values i = n,n— 1,...,1 
and all values m — 0,1, ... ,kn- 



In order to calculate P(fc„) it is useful to make the 

following definitions: 2.3 Numerical Difficulties 



/l 
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^xie ^■■ 


/2 
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^X2e 3- ■ 


/3 
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gXae^i- ■" 


fn-1 


:^ 


gairi-ie" 
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e^'. 


Equation El can now 


be written as: 


P{kn) = 


e 


'^%j['^\y) 



(13) 



(14) 



!y=0 



p(fe„ 



where f-l " is the fc„i/i derivative of /i with respect 
to y. With the above definitions, the first deriva- 
tives of the functions fi are given by: 













even milhons, 
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While the previous section gives a valid algorithm 
on how to calculate P(fc„) using equation 1141 and 
the recursive formula 1171 it turns out that the fi- 
nite precision of a normal computer will only allow 
calculations to be performed for rather small val- 
ues of fc„ before some numbers become either too 
large or too small to be stored straightforwardly in 
the computer memory. This problem is addressed 
in the following discussion. 



The factor -^j For any reasonably large number 
of dynodes, where the mean number of electrons 
coming off the last dynode, and therefore the inter- 
esting values for fc„, is typically in the thousands or 
e ^ -t:^ quickly becomes very small. 



grows to extremely large values. 



y=o 



f = 

J 71 



In- 



(.15) 

This gives a recursive formula for the first deriva- 
tive of ff. 



In order to calculate P{kn) for such values of fc„, it 

A*"" 

is necessary to absorb the small factor -^rr into the 

/j . This can be done by replacing y in equation 
1141 with py and introducing a compensating factor 



k„ 



In ~ /" ' 



(16) 



A,^ n 



P(fc„) = e--^^(- 



dy 



rfdpy) 



(18) 



■y=o 



Choosing p such that p " = -^ changes equation 
dto 



with p''" = ^ 



\fc„ 

k ' 



y=o 



(19) 



with 



Pfcn 



k„ 



HK^iy.f^'--'^ 



Pk„-lJ kn 

(fe„ 



(23) 



so the values for /^^ "^ need to be stored only for 
one iteration. 



Defining 



c "1 d 

\ 

with pfc„ 



(20) 



y=o 



ikjy 



The binomial factor When calculating /^^ " , 
using the recursive formula[52 the factor ( "7^) in 

■>kr,,i ~ 2-^ \ ^ jJk„,i-^^Jk„,i+l 



3=0 



J 



gives 



P{kn) = e-'^fiy. 



(21) 



The recursive formula established for calculating 
f{ " remains essentially the same for /^ {" : 



Jk„, 



(fc„) 



Ep"7')/:!^i-j; 



3=0 






(22) 



can get very large for large values of fc„, while 

the corresponding values for f^. i^ifk i+i §^* 
very small. To avoid the associated numerical 
problems, one can define the arrays u): j and 

(i) 
v): ^ that 'absorb' the binomial factor, such that 

equation [23 becomes: 



with /:(™j 



3=0 



(24) 



Pfe! and pfc„ 



(fcn!)' 



where 



with one additional complication. In the orig- 
inal algorithm, when calculating /j^" using the 
recursive formula El all values for /™ with 
m < kn calculated in the previous iterations^ 
could be used in the recursive formula for the 
current iteration. Now, for calculating /^^ '^' all 

values /fc ™ with m < kn, i < n have to be 
re-calculated at each iteration, because at each 
iteration the value for p in equation 1221 changes. 
To calculate /^ "" , from eauation l22l the values for 



are needed. These can be calculated using only the 
values for fl_{ j which have been calculated one 
iteration earlier: 



,U) 



^n 1 \ f*(j) 
■'k„A 



^-^^vCt-O^^- 



(25) 



2.4 Combining Results 

At each iteration A;„, before calculating /^^ ^' using 

equation 12^ the values for u^^' ■ and Vj^\, j < kn, 
are calculated from their values in the previous it- 
eration: 



-^1. 



Pk„-1 



I fcn-1 
fcn - 1 - j 



^fc„-l.' 



Jk„,i ~ Jk„-l,i 



Pk„ 



Pk„-l 






,0") 



^where P (0) , . . . , P {k„ — 1) were calculated 



j < kn 



(26) 



These results are then used to calculate: 



fe„-i 



H^) = V ^l^'.x..^-^) 



and 



.(k„ 



(fc„) ^ „*(fe„) 



starting from 



(fc„) ^ (fc„) ^ .*(fc„) ^ A^ 






and 



,(0) _ ,,(0) _ .*(0) _ f 
•'0,1 ~ "Q.i ~ Ja.i — Jt- 



(27) 

(28) 

(29) 
(30) 



where the fi are defined by equation El 



2.5 The Complete Numerical Recipe 

Using the above formulae, the problem of calculat- 
ing the probability distribution of finding fc„ elec- 
trons at the end of a PMT with n— 1 dynodes can 
be solved by a computer. A FORTRAN implemen- 
tation is listed in appendix 1X1 The program takes 
as its input the array A[n], with dimension n, which 
contains the average number of photo electrons ar- 
riving at the first dynode A[l] and the gain at each 
of the n — 1 dynodes, A[2], . . . , A[n]. The program 
fills the array ^[max] with the probabilities P[k] 
to find k electrons at the end of the dynode chain 
for all values k < max. The parameter max is also 
passed to the program. 

The values for uj^l, v)^\ needed in the recursive for- 
mulae, are stored in two two-dimensional arrays, 
where one dimension is taken by the index i ~ 
1, . . . , n, and the other by the index j = 0, . . . , max. 
As the values for ujf ^, v)^\ are needed only for one 
value of fc at a time, the arrays do not need to be 
three-dimensional; the values for u^^i,vj:- needed 
at the iteration calculating P[k] replace those from 



the previous iteration, uj^ 



(i) 



.,0-) 



The steps to calculate P[k], fc = 0, . . . , max are: 

1 Initialise program, test whether input is sen- 
sible, for example if the overall gain is larger 

than 0. Calculate all values for [^-^—) ij < 

max and store them in an array Pfracb], j = 
1, . . . max for later use. 



2 Start with calculating the probability to find 
zero electrons: k = 



3 Calculate u, 



(0) 



,(0) 



fi for i 



n, n 



!,...,!, as defined by equation 1 1 31 



- .-Ai„(0) 



4 Store the result in the array: P[0] = e 

5 Increment A; by 1. If fc > max, stop program. 



6 Calculate 



,(fc) 

'■k,n 



(k) 



hi 
k\ 



,U) JJ) 



7 Calculate MjT ^, tij, ^ for j < k and i ~ n, . . . ,1 

from u]^lii,Vf^\^ according to equation [^ 
using the values of pf|.ac[fc] calculated in step 1. 



8 Calculate u\^ ] 



„(fe) 



w^ j' for all values of i < n 
using the recursive formula 1271 Let the outer 
loop go over the index i, starting with i = 
71—1 and decrementing it by 1 until z = 1, 
and the inner loop over the summation index 
j, starting with j = and incrementing j by 1 
until j — k — 1. 



9 Store result: P[k] — e 
10 Goto step 5 



- p-Ai„(fc) 



^fc,l 



3 Fitting PMT Spectra 

3.1 Increasing Speed by Approxi- 
mating P(kn) 

When fitting PMT~pulse-height spectra, speed 
is a major problem. The number of operations 
needed to calculate P(fc„) using the recursive 
formula in eauation ll7l is 



kn i 

Nsteps ~ X! X! "-^ 
i=0 j=0 



(31) 



which becomes prohibitive for a typical PMT with 
a gain of ~ 10^ and higher. Therefore, for fit- 
ting the spectra, only the exact distribution af- 
ter the first m dynodes is calculated and then 
scaled by the gain of the remaining dynodes, ^left = 
(.9m+i.9m+2 • ■ ■ 9n-i)- When scaling the output of 
the exact distribution calculated for the first m 



dynodes, Pexact(fcm+i), to the final distribution, the 
result is convoluted with a Gaussian of width iJscaie, 
taking into to account the additional spread in the 
distribution at each remaining dynode: 



O'scah 



= \jkn 



-iCTO 



(32) 



with: 



o"o 



{gm + l9m + 2 ■ ■ ■ g-n-l) 






+ 



Sm+lS7Ti + 2 



+ ■ 



flm+l-"9n-l 



So the approximated function, P^(fc„) is 



(33) 



P^{kn)=J2- 



(j-Sleft-fcn 



J=0 



/2^^^/Jao 



KV?-o) p(^^y (34) 



In practice the sum only needs to be calculated for 
values of j ■ g\eft that are a few Uscaie around A:„. 



3.2 Reducing the Number of Param- 
eters 

P(fc„) depends on n parameters: the gain of each 
dynode and the number of photoelectrons produced 
in the cathode. For the case of the 12-dynode 
PMT, there are 13 parameters. It is possible, how- 
ever, to reduce this number to two: 



3.3 Adapting the Function to Fit 
Measured Data 



In practice, spectra are not measured in numbers 
of photoelectrons, but in ADC counts digitised by 
the readout electronics. The function describing 
the spectra needs to relate the ADC counts, fcadc, 
to the number of electrons at the end of the dynode 
chain, fc„. This requires two parameters: the offset, 
or pedestal mean, poi and the conversion factor, c„ 
of kn to ADC counts. The resulting function is 
convoluted with a Gaussian of width a to take into 
account electronics noise: 



-f*cont(ri^adc) 



27rcr 



* (-P ((fcadc - Po) /Cn) ■ C„) , 



(36) 



where * is the convolution operator. -Fcont treats 
fcadc as a continuous variable, with a one-to-one 
relation between /cadc and fc„; in fact the readout 
electronics deliver only integer-value ADC counts, 
integrating over the corresponding pulse heights. 
Thus the final function for describing ADC spectra 
is: 



/■Kadc + ll-O 

^(fcadc) = / ^cont(fc^dc)*:dc 

"'fcadc-0.5 



(37) 



1. the mean number of photoelectrons produced 
in the photo cathode 



4 Example Fits 



2. the gain at the first dynode. 



Using 



gcxF", 



(35) 



where V is the voltage difference over which the 
electron is accelerated, the gain at the other dyn- 
odes can be calculated from the gain at the first 
dynode. The parameter a has values typically 
between 0.7 and 0.8 |HamOO) : in the following, 
a — 0.75 is used. 



The fits are performed as binned log-likelihood fits: 
for each 1-ADC-count wide bin fcadc, containing rii 
events, the binomial probability of having rii "suc- 
cesses" in A'aii trials is calculated, where A'aii is the 
total number of events. The probability of an indi- 
vidual "success" is given by P ( fcadc ) • 

The probability distribution for the number of elec- 
trons after the fourth dynode is calculated without 
approximation. Then the function is scaled, ap- 
proximating the additional spread due to the re- 
maining dynodes with a Gaussian, as described in 
the previous section. 



Table 1: Voltage distribution in 12-dynode PMT, normalised to the voltage between dynodes 3 and 4. 



voltage 
dynode number 



|3|2|2|1|1|1 
Cathode I 1 I 2 I 3 I 4 I 5 I 



1 I 1 I 2 I 
1 10 111 1 12 



Figm-e 1: MC-generated PMT ADC-spectrum, Table 2: Monte Carlo input compared with mean 



from 100k events, with Ai 
perimposed 



0.15. The fit is su- and RMS of the results from fits to 128 simulated 
spectra, with Ai = 0.15 



I 
u 

Q 
< 



10 



10 



10 



10 



10 





■ 
































^ 




^. 


L 










V/^ 


s 


=5* 




bj 


!^ 






^ 



75 100 125 150 175 200 225 250 275 

ADC counts 



4.1 MC-Generated Spectra 

The validity of the the method has first been es- 
tablished on Monte Carlo simulated data. The 
Monte Carlo program simulates the output of a 
PMT pixel. The gain at the first dynode is gi = 5 
and the gains at the other dynodes are calculated 
from g oc V^ with a = 0.75. The values for V are 
given in tabled 

The fit function is applied to two sets of 128 simu- 
lations with 10^ events each, one set with 0.15 pho- 
toelectrons per event, one with 3.0 photoelectrons 
per event. A spectrum from each set is shown in 
figures n and |21 

The fits are performed varying the gain of only one 
dynode and calculating the gains at the other dyn- 
odes using the same value for a as in the Monte 
Carlo program that generated the spectrum. The 
fit results agree very well with the input values, as 
shown in tables El and 13 To test the sensitivity of 
the fit result on the exact knowledge of a, the fit to 
the spectrum in figure ^ is repeated assuming dif- 





MC input 


Mean fit result 
± RMS spread 


Ai 


0.150 


0.1501±0.0013 


.91 


5.000 


5.0012±0.058 


Po 


100.00 


99.999±0.0038 


a 


1.0000 


1.0004±0.0027 


c« 


3.20- 10""* 


(3.23±0.46) • 10-4 



Figure 2: MC-generated PMT ADC-spectrum, 
from 100k events, with Ai = 3. The fit is super- 
imposed. 
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Table 3: Monte Carlo input compared with mean 
and RMS of the results from fits to 128 simulated 
spectra (representing 2 64-channel MaPMT's), 
with Ai = 3 





MC input 


Mean fit result 
± RMS spread 


Ai 


3.000 


3.002±0.022 


.91 


5.000 


4.985±0.107 


Po 


100.000 


99.999±0.021 


a 


1.000 


0.999±0.016 


Cn 


6.4- 10--* 


(6.45±0.17) • 10--* 



Table 4: Monte Carlo input compared with fit-result for the MC-generated pulse height spectrum shown 
in figure ^ using different assumptions in the fit. 





MC input 


Fit result 


Fit result 


Fit result 


I'it result: 3 




Q = 0.75 


Q = 0.75 


Of = 0.5 


a = 1 


indep. dyn's 


Ai 


0.1500 


0.1490 


0.1491 


0.1489 


0.1492±0.0013 


31 


5.00 


5.039 


4.852 


5.291 


4.74±0.44 


52,53,512 


51 • (f )'^ = 3.69 
gi •(!)" = 2.19 








4.51±1.35 


54, ■ • ■ ,511 








1.97±0.21 


Po 


100.000 


100.000 


100.000 


100.000 


100.000±0.003 


a 


1.0000 


1.0028 


1.0029 


1.0027 


1.0028±0.0025 


c„ 


3.20 ■ 10"^ 


2.90 ■ 10"'' 


0.373 ■ 10-* 


19.7- 10-* 


(4.37 ±0.98) • 10-* 



ferent values for this parameter in the fit-function: 
a = 0.5 and a = 1.0. The results are given in table 
0] Another fit was performed that does not use the 
formula g oc V"' . Here it is only assumed that dyn- 
odes with the same accelerating voltage have the 
same gain. Instead of one gain, three gains need 
to be fitted, one for each accelerating voltage. The 
fits are performed using the function minimisation 
and error analysis package MINUIT |.Tam94| . The 
results from this fit, with error-estimates provided 
by MINUIT, are given in the last column of table 

H 



Figure 3: Data from 6k events in black, with fit 
superimposed. The dashed line indicates the sin- 
gle photoelectron contribution. The signal loss 
refers to the fraction of photoelectrons lost below 
the threshold cut; both the total fraction of photo- 
electrons, and the fraction of single photoelectron 
events lost below the cut is given. These numbers 
do not include the loss at the first dynode due to 
photoelectrons producing zero secondary electrons. 



Comparing the results for the different assumptions 
shows that they have little impact on the the fit- 
ted value for the number of photo electrons and the 
gain at the first dynode. Most of the error intro- 
duced by an incorrect estimate of the parameter a 
is absorbed into the ratio of ADC~counts to elec- 
trons, c„, while the values for Ai and 171 come out 
close to the input values. 



4.2 Application to Testbeam Data 
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The fit method has been applied to spectra ob- 
tained from a prototype RICH detector, incorpo- 
rating an array of nine 64-channel Hamamatsu 
PMTs and operated in a CERN testbeam RICOOJ . 
Fits were performed to estimate the signal loss at 
the first dynode and below the threshold cut. 

Figure |3| shows an example of such a fit to a spec- 
trum obtained in the testbeam. The fit describes 



Table 5: Result of fit applied to testbeam data 





Fit result 


Ai 


0.107±0.005 


.91 


3.60±0.20 


Po 


43.06±0.01 


a 


0.724±0.008 


c„ 


(61±30) • 10-'' 



Figure 4: MC-simulated PMT spectrum with back- 
ground; the fit result is superimposed 
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MC input 


Fit result 




^^1 
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Figure 5: The different contributions to the fit in 
figure 
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4.3 Background 

Apart from the Gaussian noise taken into account 
here, various other sources of background, such as 
electrons released due to the photoelectric effect in 
the first dynode, thermal electrons from the pho- 
tocathode or the dynodes, genuine photoelcctrons 
missing the first dynode, etc, can contribute to a 
PMT pulse height spectrum. A detailed discus- 
sion of such background is beyond the scope of this 
paper. However, any type of background that orig- 
inates from within the dynode chain can be natu- 
rally accomodated in the fit method described here, 
since this background undergoes the same type of 
amplification process as the signal. To illustrate 
this, a spectrum has been generated with a Monte 
Carlo program, assuming a signal of 1.6 photoelcc- 
trons per event and a background of 0.16 photo- 
electrons per event due t o the ph otoelectric effect 
in the first dynode (see [CZ+Oll for a fit to real 
data showing this kind of background, using a dif- 
ferent method) . The function to fit this spectrum is 
obtained by convoluting the background-free func- 
tion P(A:) with another function Pbg(fc). Pbg(fc) is 
identical to P(fc) except that the amplification due 
to the first dynode is missing and that the num- 
ber of photoelcctrons per event hitting the second 
dynode, Abg, is a new free parameter. In the exam- 
ple given here, P and Phg{k) are calculated to give 
the exact distributions for signal and background 
respectively after the fourth dynode; then the two 
distributions are convoluted and the result is scaled 
according to eauation l34l The generated spectrum 
and the fit result are shown in figure^] the fit func- 
tion is shown again in figure El showing the non- 
pedestal and the single photoelectron contributions 
separately. 



5 Summary 



the data well, with a x^/dgf of 1.2^. The line in 
figure 121 marks the threshold cut used for photon 
counting in the testbeam. The fraction of single 
photoelectron events below that cut is ~ 10% (this 
does not include the irrecoverable loss of photoelcc- 
trons that do not produce any secondaries in the 
first dynode). 



An analytical formula for the the probability dis- 
tribution of the number of electrons at the end of 
a dynode chain, or any "snowball" like process de- 
scribed by a series of Poisson distributions, is de- 



^The fit is performed with the same log— UkeUhood 



method that was used for the MC spectra; 
calculated after the fit. 



a X value is 



rived. The formula describes the amplification pro- 
cess at all stages exactly, in particular without ap- 
proximating Poisson distributions with Gaussians. 
ft is evaluated as a function of the number of 
photoelectrons coming from the cathode and the 
gains at each dynode. The initially found formula 
is adapted to reduce numerical problems due to 
the multiplication of very large numbers with very 
small ones. A numerical recipe is given that imple- 
ments that function. 
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FORTRAN Routine to Cal- 
culate P{kr, 



It is shown how the function can be used as the core 
element of an approximated, but faster algorithm, 
that calculates the exact distribution for the first 
few dynodes and then scales the result according 
to the gain at the remaining dynodes, approximat- 
ing the additional spread at those dynodes with a 
Gaussian. The number of dynodes for which the 
distribution is calculated exactly is not limited in 
principle and can be adjusted according to the pre- 
cision required, and the computing time available. 
It is also shown how to modify the function to de- 
scribe ADC-spectra obtained from read-out elec- 
tronics, rather than directly the number of elec- 
trons at the end of a dynode chain. 



This fast algorithm is then used to fit Monte Carlo 
generated ADC-spectra. In the fit function, the 
electron distribution after the first four out of 
twelve dynodes is calculated exactly. The fit re- 
sults reproduce the MC-input values well. The de- 
pendence of the fit result on the assumptions made 
to reduce the number of fit-parameters is inves- 
tigated. These results show that the fitted value 
for the number of photoelectrons per event is very 
weakly dependent on the different assumptions con- 
sidered here, and the fitted gain on the first dyn- 
ode also does not depend strongly on them. Real 
data from a multi-anode PMT used in the 1999 
LHCb- RICH testbeam are fitted, and shown to be 
described well by the function. Finally it is illus- 
trated how the fit function can be modified further 
to accommodate background from within the dyn- 
ode chain, using the example of the photoelectric 
effect in the first dynode. 



SUBROUTINE DYNODE_CHAIN(OUT, MAX, LAMBDA 

IMPLICIT NONE 

This program takes as its input the maxii 

the end of the dynode chain, for which i" 

MAX , the average number of phi 

dynode, LAMBDA(l), the gains i 

... LAMBDA (DYNODES) and the d 



DYNODES) 



m number of electrons at 
should calculate PCk_n), 
hitting the first 
lach dynode, LAMBDA(2), 
ision of the array LAMBDA: 



DYNODES. It calls the routii 
end of this file. 



I MAKE_P_RATIO, which is listed at the 



The output is put into the array OUT(MAX) , where the probability 
to find k_n < MAX electrons at the end of the dynode chain is 
given by OUT(k_n). 

Written by Jonas Rademacker. 

INTEGER MAX, DYNODES 

DOUBLE PRECISION OUT(0:MAX), LAMBDA (DYNODES) 

INTEGER ABS.MAX, MAX.DYN 
PARAMETER(ABS_MAX=50001,MAX_DYN=13) 

INTEGER IX,IY,M,I, K, J 

To avoid having to define a limit on the number k_n that can be 
calculated, one could create these arrays outside the program and 
pass them on. 

DOUBLE PRECISION F(1:MAX_DYN) ! corrsponds to f-{star> in the text 
DOUBLE PRECISION U(1:MAX_DYN,0: ABS.MAX) ,V(1:MAX_DYN,0: ABS.MAX) 
DOUBLE PRECISION X(i:MAX_DYN) 

DOUBLE PRECISION FASTNULL 
PARAMETER (FASTWULL=l.d-300) 



DOUBLE PRECISION MEAN 

DOUBLE PRECISION P_rat 

INTEGER MAX.OLD 

SAVE MAX. OLD 

DATA MAX. OLD/ -9999/ 

SAVE P.ratio 



d(ABS_MAX), F.FACTOR, U.FACTOR, V.FACTOR 



— Some initialisations and tt 
DO IX=1,MIN(ABS.MAX,MAX) ,+1 

OUT(IX)=0.dO 
ENDDO 
IF (ABS.MAX . LT . MAX) THEN 

RETURN 
ENDIF 

MEAN = l.DO 
DO IX=1, DYNODES, +1 

MEAN = MEAN*LAMBDA(IX) 
ENDDO 
IF(MEAN.LE.O.dO)THEN 

OUT(0)=l.dO 

RETURN 
ENDIF 

— make and save the factors I 
IF(MAX.GT.MAX.OLD)THEN 

MAX.OLD=MAX 

CALL MAKE_P_RATIO(P_ratio,^ 
ENDIF 

— Calculate the probability 1 
F(DYNODES)=l.dO 



.io(k) = (p.{k>/p.{k-l»-tk> - 



) electrons (k_i 
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U(DYNDDES , 0) =FCDYNDDES) 
V (DYNDDES , 0) =F (DYNDDES) 
DO IX=DYNDDES-1,1,-1 

X(IX) = LAMBDA(IX)*DEXP(-LAMBDA(IX+1)) 

F(IX) = DEXP(X(IX)*F(IX+1)) 

U(IX,0) = F(IX) 

V(IX,0) = F(IX) 
ENDDD 
DUT(0)=DEXP(-LAMBDA(1))*F(1) ! < save the result 

— Calculate the probabilities for k_n=l MAX elei 

DO K=1,MX,+1 

calculate f_ii 

IF(F(DYHDDES) . LT . FASTNULL) THEN 
F(DYNODES)=0.dO 

ELSE 

F(DYMODES)=F(DYNODES) * LAMBDA (DYNDDES)/DBLECK) 

ENDIF 

U(DYNDDES,K)=F(DYHDDES) 

V(DYNDDES,K)=F(DYHDDES) 

re-calculate U and V from previous iteration; 
DO J=0,K-1,+1 

F_FACTOR=P_ratio (K) ** (DBLE( J) /DBLE(K) ) 
IF(K-1-J.GT.0)THEN 

U_FACTOR=DSqRT(DBLE(K-l)/DBLE(K--i-J))* 
& F.FACTOR 

ELSE 

U_FACTOR=F_FACTOR 
ENDIF 

V_FACTOR=DSqRT((DBLE(K-l)/DBLE(K-J)))* 
& F.FACTOR 

DO 1=DYN0DES,1,-1 

U(1,J)=U(1,J)*U_FACT0R 

V(I,J)=V(1,J)*V_FACT0R 
ENDDO 
ENDDD 
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apply the recursive formula to get f~tk}-_i 
DO I=DYN0DES-1, 1, -1 
F(l)=0.dO 
DO J=0,K-1 

F(I)=F(1)+U(1,K-1-J)*X(1)*V(1+1,J+1) 
ENDDO 

UCI,K)=F(1) 
V(I,K)=F(1) 
ENDDO 
. calculate P(k): 

0UT(K)=DEXP(-LAMBDA(1))*F(1) ! < save the ] 

ENDDD 

RETURN 
END 



SUBROUTINE MAKE_P_RAT10(P_ratio,MAX) 

IMPLICIT NOME 

INTEGER MAX 

DOUBLE PRECISION P_ratio(MAX) 

INTEGER N 

DOUBLE PRECISION NFAC 

DOUBLE PRECISION PI, E 

PARAMETER (PI =3 . 1415927d0 , E=2 . 718281828d0) 

INTEGER APPROX.FROM 
PARAMETER (APPR0X_FR0M=25 ) 

NFAC=1.D0 

P_ratioCl)=l.dO 

DO N=2,MIN(APPR0X_FR0M-1,MAX) ,+1 

NFAC=NFAC*DBLE (N- 1 ) 

P_ratio(N)=(NFAC**(l.dO/DBLE(N-l)))/DBLE(N) 
ENDDO 

DO N=APPR0X_FR0M,MAX,+1 
P_ratio(N)= 
& (2.D0*PI*DBLE(N-1))**(1.D0/(2.D0*DBLE(W-1)))* 
& DBLE(N-1)/(E*DBLE(N))* 
& (I.d0+I.d0/DBLE(12*(N-1))+ 
& I.d0/DBLE(288*(N-1)**2) 
k )**(1.D0/DBLE(N-1)) 

ENDDO 

RETURN 
END 
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